The immune cell atlas of human neuroblastoma

Summary Understanding the complete immune cell composition of human neuroblastoma (NB) is crucial for the development of immunotherapeutics. Here, we perform single-cell RNA sequencing (scRNA-seq) on 19 human NB samples coupled with multiplex immunohistochemistry, survival analysis, and comparison with normal fetal adrenal gland data. We provide a comprehensive immune cell landscape and characterize cell-state changes from normal tissue to NB. Our analysis reveals 27 immune cell subtypes, including distinct subpopulations of myeloid, NK, B, and T cells. Several different cell types demonstrate a survival benefit. In contrast to adult cancers and previous NB studies, we show an increase in inflammatory monocyte cell state when contrasting normal and tumor tissue, while no differences in cytotoxicity and exhaustion score for T cells, nor in Treg activity, are observed. Our receptor-ligand interaction analysis reveals a highly complex interactive network of the NB microenvironment from which we highlight several interactions that we suggest for future therapeutic studies.

INTRODUCTION Neuroblastoma (NB) is a pediatric cancer deriving from the neural crest during embryonic development and arises primarily in the sympathetic nervous system. 1 It has a wide range of clinical presentations and outcomes, from spontaneously regressing tumors to highly aggressive and metastatic cases. This heterogeneity is reflected by large differences in 5-year survival rates ranging from 90% to 95% for low-and intermediate-risk disease, to 40%-50% for high-risk disease. 2 The role of the immune system is widely recognized as critical in cancer development, progression, and therapy resistance in adults. Various immune cell subsets can either directly or indirectly support or suppress tumor growth in adult cancer. 3,4 Multiple studies in NB have pointed to a critical role for the immune system in both prognosis and response to treatment. For instance, the infiltration of T cells in therapy-resistant NB is known to improve clinical outcome, 5 while the prognostic value of the specific T cell subtypes remains unclear. 6 Ex-vivoexpanded tumor-infiltrating lymphocytes from patients showed phenotypic heterogeneity; however, they were non-reactive toward autologous tumor cells. 7 In addition, T cell receptor sequencing on NB-derived T cells revealed clonal expansion in only a small number of untreated patients, [8][9][10] indicating that a tumor antigen response may occur only in a limited number of patients. Furthermore, T cell memory formation, and a sug-gested active tumor microenvironment, have also been described in NB, 11 implying an active inflammatory response in these tumors. However, the recently described tissue-resident T cells, which highly resemble T memory cells, demonstrated that a memory T phenotype does not have to indicate tumor antigen recognition and memory formation. 12 Other lymphocytes, such as invariant natural killer (NK) T cells have been implied to play an important role in the spontaneous regression of NB due to their prominent presence in low-risk tumors. 13 Also, NK cell marker gene and protein expression are correlated with improved prognosis, 14 and infiltration of NK cells into NB tumors has been detected by bulk RNA sequencing. 15 Furthermore, different myeloid cell subtypes play an important role in disease progression and metastasis in NB, 16,17 where they are described as highly heterogeneous. 18,19 An increase in the number of macrophages was detected in metastatic NB patient samples compared with local primary tumors, 17 and macrophages are associated with poor prognosis 20 ; whereas dendritic cell gene and protein programs in NB tumors in turn have been correlated with improved prognosis. 14 In addition, the presence of myeloidderived suppressor cells has been observed in both human and mouse NB tumors and has been associated with poor prognosis. 21,22 These insights into the function of immunity in NB have led to the development of anti-GD2 immunotherapy that is given in conjunction with conventional treatment in high-risk NBs. 23  Although immune cell populations have been studied, many of the studies have either conveyed conflicting results, focused on a single immune cell type, or used a limited number of patient samples. 6 Previous single-cell RNA sequencing (scRNA-seq) studies in NB have mainly focused on the tumor cell compartments [25][26][27][28][29][30][31] and the cells in the tumor microenvironment form a complex interactive network of cells and molecules, where the immune cell function is dependent on the different cells working together in the tissue. Therefore, providing a systems biology characterization of the immune cells, their prognostic impact, and their interactions with tumor cells, and each other, at the single-cell level, will provide an improved understanding of the NB microenvironment that can be exploited for future therapies. In this study, we performed scRNA-seq on human NBs, we contrasted our data with a single-cell dataset from normal fetal adrenal tissue and additional NB datasets to provide an in-depth analysis and characterization of the different immune cells, their molecular state shifts, and cellular networks.

Global immune cell landscape of human NB
We performed single-cell transcriptomics (10X Chromium Single Cell Solution, see STAR Methods) on tumor tissue from 17 NB patients (19 samples). Tumor specimens were dissociated using a protocol that enriches for immune cells ( Figure 1A). 32 The samples were derived from all clinical risk groups and included patients with both primary and recurrent disease (Table S1). A preprint article by our lab focused on the tumor and stroma compartment, 28 whereas in this study we focused on the immune cells of NB. We used the algorithm CONOS 33 for joint analysis of all NB samples and created an immune cell landscape ( Figure 1B; for cellular map including tumor and stroma cells, see Figure S1A). Our initial analysis of the immune cells including 46,261 cells revealed a rich repertoire of the main immune subtypes, including myeloid, B, T, and NK cell lineages, and consisted of 10 clusters in total ( Figures 1B-1D). The global immune cell clusters were identified using key marker genes: myeloid cells (LYZ, C1QA, S100A9, and CD1C), 34,35 B cells (MS4A1, BANK1, and CD79A), 36,37 T cells (CD3D), 38 and NK cells (CMC1 and GNLY) 39 (Figures 1C and 1D). We validated immune cell infiltration of the different immune cell types detected in our scRNA-seq data by performing immunohistochemistry on a separate cohort of 43 patients (Table S2) taking one or multiple samples from the same patient ( Figure S1B). We stained for key markers detecting T cells (CD3+), B cells (CD19+), NK cells (NKp46+), macrophages (CD68 + ), dendritic cells (CD1a+), and neutrophils (NE+) ( Figure 1E). Neutrophils were absent in the scRNA-seq data, reflecting a known limitation of the single-cell method utilized. 40 The majority of immune cells detected in the single-cell analysis were T cells and myeloid cells where proportions for each immune cell subtype differed between patients (Figures 1F, S1C, and S1D). In addition, we detected T cells and macrophages for most of the samples we performed immu-nohistochemistry on, and immune cell infiltration overall was between 5% and 35% of total cells in the tumor (Figures 1G and S1E). Next, we created a combined dataset containing data from previous NB scRNA-seq studies [25][26][27] and we determined the fraction of immune cells in low-, intermediate-, and highrisk disease. We detected no differences in the major immune cell populations comparing the risk groups ( Figure S1F). Survival analysis on bulk RNA sequencing data for 498 samples, based on gene expression signatures derived from our scRNA-seq (Table S3), showed that high expression of immune-related genes was significantly correlated with improved survival in human NB ( Figure 1H).

Infiltration of myeloid cells with distinct cell states detected in NB
Myeloid cell infiltration has been described for multiple cancers and shown to support tumor growth. 18,41,42 In NB, discrepancies between studies and a lack of study into myeloid cell heterogeneity in NB prompted us to analyze these cells in detail.
We have used in silico subcluster analysis 32 focusing on the myeloid cells present within NB tissues obtained by our lab to identify nine distinct subpopulations. These we annotated as Mono-1/2, CLEC9A+ myeloid dendritic cells (mDCs), CD1C+ mDCs, mature-LAMP3+ mDCs, and Macro-1, 2, 3, and 4, based on their expression of key marker genes for respective cell lineages (Figures 2A-2C). Multiplex immunohistochemistry, on patient samples from the single-cell cohort (Table S1), was used for the detection of antigen-presenting myeloid cells (CD11c+, five out of five patients, and HLA-DR+, in three of three patients) ( Figures S2A and S2B). To characterize the cell state of the myeloid cells, we curated a gene signature score based on existing scRNA-seq studies describing previously characterized tumorderived monocyte and macrophage cell states (Table S4). Our analysis revealed a significantly higher monocyte score in both Mono-1 and Mono-2 compared with the other myeloid populations ( Figures S2C and S2D), substantiating that these cells are monocytes. Macro-2 and Macro-3 also had a high monocyte score ( Figures S2C and S2D). Furthermore, Macro-1 showed the highest macrophage cell identity score followed by the three other macrophage populations when compared with the monocytes (Figures S2E and S2F), substantiating their macrophage cell state.
Focusing on Mono-1/2 subclusters, we annotated these clusters based on the expression of the key monocyte marker genes CLEC7A, NLRP3, and BST1. 43 Both monocyte clusters expressed the pro-inflammatory cytokine IL1B, the monocyte genes TIMP1, CD44, and G0S2, and several inflammatory-related S100A genes ( Figures 2B and 2C). Mono-1 appeared to consist of two subpopulations within: a small population of non-classical monocytes with expression of FCGR3A, FCGR3B, and IFITM3, and a larger population of cells that exhibited specific expression of TIMP1, S100A4, S100A12, and CD55 ( Figures 2B and 2C). Mono-2 had a mixed gene expression pattern with similarities to both monocytes and macrophages. This population might be differentiating into Macro-2 since part of the population displayed expression of the  2C). 35 To provide a possible function of Mono-1/2, we applied a functional gene signature score analysis focusing on the inflammatory activity of monocytes, and compared them to macrophages (Table S4). Our analysis showed that Mono-1/2 have a significantly higher inflammatory monocyte score compared with macrophage populations and are therefore most likely pro-inflammatory monocytes (i.e., TIMs) ( Figures 2D and 2E) that we previously characterized in human prostate cancer bone metastases. 32 Further investigation of the macrophage populations showed that the macrophage clusters had high expression of the classic macrophage marker gene CD68 ( Figure 2B). 35 Furthermore, Macro-1 displayed a population-specific gene expression pattern of tumor-associated macrophage markers, such as TREM2, FABP5, and FABP4 ( Figures 2B and 2C). [44][45][46] Macro-2 expressed high levels of neutrophil-attracting chemokines, including CXCL2, CXCL3, and CXCL8. 47 In addition, this population exhibited high expression of GPR183, DUSP2, CD83, and HLA class II genes ( Figures 2B and 2C). Macro-3 displayed expression of other macrophage genes, such as LYVE1, CSF1R, and MRC1, 20,45,48 and Macro-4 had the highest level of HLA class II genes ( Figures 2B and 2C). To demonstrate the possible functional properties of the different macrophage populations, we ran analysis of M1 (i.e., immune active) and M2 (i.e., immunosuppressive) signature score (Table S4), which revealed that Macro-1 had a significantly higher M2 score than the other macrophage populations, and those in turn had a higher M2 score compared with the monocytes (Figures 2F and 2G). Among Macro-1 to 4, Macro-2 had the highest M1 signature score ( Figures S2G and S2H).
We next examined the myeloid dendritic cell population that branched into three different subpopulations, including CLEC9A+ cells, CD1C+ cells, and mature-LAMP3+ cells. CLEC9A+ mDCs expressed high levels of CLEC9A, RAB7B, BATF3, WDFY4, and CADM1, which phenotypically corresponded to DC1 cells (Figures 2A-2C). 43 CD1C+ mDCs were identified by the high expression of CD1C, CLEC10A, and FCER1A, but also exhibited high expression of HLA class II antigen-presenting genes, and therefore had a similar phenotype to previously described DC2. 43 Mature-LAMP3+ mDCs showed high expression of CCR7 and CD83 ( Figures 2B and 2C).
To increase the power of our analysis, but also to compare our findings with corresponding normal tissue for NB, we combined our dataset with three recently published scRNA-seq studies (including a recent study 25 ) of normal fetal adrenal tissue (i.e., where NB is believed to originate) and NB tumors (see section Data and code availability in STAR Methods). [25][26][27] Here, we observe that the fetal adrenal data mainly contained myeloid cells ( Figures S2I and S2J). In the combined dataset, we detected the same myeloid subpopulations as in our dataset, validating our findings ( Figures 2H and S2K), where the annotations were based on key marker genes ( Figure S2L). To determine the role of specific myeloid populations on disease risk stratification, we compared the fraction of cells and the different functional gene signature scores from the combined dataset [25][26][27] in low-, intermediate-, and high-risk NB compared with healthy fetal adrenal myeloid cells. Our analysis revealed that the main myeloid population found in fetal adrenal tissue is Macro-3. Also, these cells showed significantly lower cell fractions in NB patients, while Macro-4 demonstrated an increase in the intermediate-and high-risk NB ( Figures 2I and 2J). In addition, Macro-3/4 were close together in the embedding (Figure 2H). Looking into cell state, the different gene signature scores showed a significantly increased M1 score for four out of six populations, and a significantly decreased M2 score with higher disease risk in Macro-3 and Mono-2 ( Figures S2M and  S2N). Surprisingly, we discovered a significant increase of inflammatory monocyte score in all populations except Macro-4, indicating that the presence of tumor cells strongly initiates inflammatory gene expression in these cells ( Figure 2K). Focusing on Mono-1 and Mono-2 in this scoring analysis, we also detected a difference in gene expression comparing fetal adrenal and NB cells ( Figure 2L).
Finally, to evaluate for potential prognostic value of different myeloid cell populations, we ran survival analysis on bulk RNA sequencing data based on key marker gene expression from our scRNA-seq dataset (Figures 2C; Table S5). Interestingly, Mono-2 and Macro-4 significantly correlated with improved survival, whereas Macro-1 was significantly correlated with decreased survival ( Figure 2M). Signatures based on CLEC9A+ and CD1C+ dendritic cell populations were associated with improved survival, whereas the other monocyte and macrophage populations and mature mDCs did not show an association with survival ( Figure S2O). Since NB is a clinically heterogeneous disease, we split the bulk RNA sequencing data to determine survival into low-and high-risk cases, and into non-MYCN-amplified and MYCN-amplified cases. Next, we applied the gene scores for the different myeloid cell populations on the split data and created a summary figure. There, we detected a significant correlation with improved survival for CD1C mDCs in the low-risk group, but there were no differences for the other groups ( Figure S2P).

B and NK cell heterogeneity found in NB
We next focused our analysis on the B and NK cells. Subcluster analysis of the B cells in the combined dataset revealed the presence of four subpopulations in NB tissues ( Figures 3A-3C). B cell infiltration was confirmed by multiplex immunohistochemistry on patient samples from the single-cell cohort (Table S1, Figures S3A and S3B). Active B cells expressed B cell marker genes, such as IGHD, IGHM, CD69, and CD83, 49 and plasma cells expressed high levels of MZB1 and JCHAIN (typical B cell plasma marker genes) ( Figures 3B and 3C). 50,51 In addition, we detected memory B cells that were annotated based on the expression of CD27, CLECL1, and ZBTB32 ( Figures 3B and 3C). 52 We also   Figures 3D and S3C). In addition, we determined the prognostic value of the different B cell subtypes in a similar manner as performed in Figure 1H and for the myeloid subpopulations (genes can be found in Figure 3C and Table S5), where we observed a trend in the presence of plasma cells and tertiary lymphoid structure (TLS) signature score with survival; however, the results were not statistically significant ( Figures S3D and S3E). Finally, in a similar approach as for the myeloid cells, we created a survival summary figure showing the different clinical conditions, where we detected no significant correlation with survival ( Figure S3F).
NK cells have potent tumor-killing properties, and NK cell infiltration has been demonstrated as a good prognostic marker in different cancers. 56,57 In the combined NB dataset, subcluster analysis of the innate lymphoid cells (ILCs) identified four subclusters ( Figure 3E). No NK cells were detected in the fetal adrenal glands (data not shown). The NK-active cell cluster was annotated based on the expression of immediate-early genes, such as ZPF36, NR4A2, JUNB, FOS, FOSB, and CD69. 39 Classical NK cell subsets CD56 bright and CD56 dim were detected ( Figures 3F and 3G). 39,58 In addition, an ILC3 cell cluster was annotated based on their lack of expression of classical NK cell marker genes TBX21 and EOMES and the high expression of AHR, KIT, IL7R, and RORC ( Figures 3F and 3G). 59 The potently cytotoxic CD56 dim NK cells ( Figure S3G) were found significantly enriched in cell number in intermediate-and high-risk disease ( Figures 3H, 3I, and S3H). However, their presence in NB did not correlate with improved survival ( Figure S3I). On the other hand, active NK cells were significantly correlated with improved survival ( Figure 3J). No significant differences were detected when the data were split by risk group and MYCN status ( Figure S3J). NK cell survival analysis was done as described above using key marker gene expression (Figure 3G; Table S5).
Distinct subtypes of T cell infiltrates correlate with improved NB survival Various subpopulations of T cells have been suggested to play a critical role in cancer development and progression, which has prompted the identification of new immuno-therapeutic approaches based on modulating T cell function and dysfunction in cancer. [60][61][62][63] In our NB dataset, subcluster analysis revealed ten distinct subclusters of T cells, including two CD4 + clusters (Tregs [regulatory T cells] and Th17), four cytotoxic T lymphocyte (CTL) CD8 + clusters (CTL-1/4), a naive T cell cluster containing cells with CD4 expression and cells with CD8A/B expression, a proliferating T cell cluster, and an NKT cell cluster ( Figures 4A-4C). Infiltration of the different T cell subtypes was validated by multiplex immunohistochemistry (patient information in Table S1) and single staining immunohistochemistry (patient information in Table S2) (Figures S4A-S4D). Single-cell transcriptomics and immunohistochemistry both showed the presence of Th17 cells (RORC, CCL20, CCR6, and, less prominently, IL17A). 64 However, other cell types, such as ILC3 59 and gdT cells, 65 can also produce IL-17, where the latter promote tumor cell proliferation and migration in NB cell lines. 65 Tregs were identified based on FOXP3 and IL2RA expression, two key Treg marker genes 66 (Figures 4B and 4C), and infiltration of Tregs was validated by FOXP3 staining (Figures S4C and S4D). Because of their tumor-promoting role in adult cancer, we sought to evaluate the activity of Tregs in NB based on a gene signature score curated from literature (Table S4). Treg activity was evaluated in the combined dataset, where we detected no difference in Treg activity score comparing low-, intermediate-, and high-risk disease ( Figures S4F and S4G). Surprisingly, the number of Tregs showed a significant decrease from intermediate-to high-risk disease ( Figure 4D). In addition, Tregs did not correlate with survival ( Figure S4H). Due to the low number of T cells in the normal fetal adrenal data this was not included in the combined data analysis.
Next, we sought to delineate the CD8+ T cell compartment, which revealed four different subtypes of cytotoxic T lymphocytes where CTL-1 cells expressed high levels of PRF1 and GZMB ( Figures 4A-4C). The expression of granzyme-B and perforin-1 showed high cytotoxic properties for CTL-1. CTL-4 expressed high levels of TNF, CD69, and IFNG. Cytotoxic T cells upregulate the expression of CD69 and IFNG upon TCR stimulation, 67,68 showing that CTL-4 are activated T cells (Figure 4B). To further investigate the functional properties of the CTL populations in NB, we took the combined dataset (Figure S4E) and performed cytotoxicity and exhaustion signature score analysis comparing low-, intermediate-, and high-risk disease (Table S4). Interestingly, no differences were found for both cytotoxicity and exhaustion score for the different CTL populations in NB (Figures S4I and S4J).
Since different T cell subtypes have rendered diverse results in prognostic value in different NB studies, 5,15,69,70 we ran survival analysis for the different T cell subtypes based on their marker genes ( Figure 4C; Table S5). Including all data, CTL-1, CTL-3, CTL-4, naive T cells, and Th17 signatures were significantly correlated with improved survival (Figure 4E), whereas proliferating T cells, CTL-2, and NKT cell signatures did not ( Figure S4K). When we split the data into risk group and divided by MYCN status as described for other cell types above, we detected no significant differences ( Figure S4L).

Ligand-receptor interaction analysis reveals a highly complex interactive network of immune cells with tumor and stroma
To investigate potential interactions and molecular crosstalk between cells located within the NB microenvironment, Cell Reports Medicine 3, 100657, June 21, 2022 7 Article ll we performed receptor-ligand interaction analysis. 71 To find interactions that may be targetable, we used significant correlation to survival as a cutoff. In addition, we set a cutoff for specific expression of both ligand and receptor in a specific cellular subtype. Our analysis revealed that the tumor microenvironment of NB consists of a highly complex interacting cellular network ( Figure 5A). Most of the interactions involving tumor cells occurred with the different stroma cell populations, but there were also noticeable interactions within tumor/stroma cells and the immune cells ( Figures S5A and S5B). Importantly, we detected a high number of potential interactions between stroma and immune cells, involving CTL1, mesenchymal, and endothelial populations ( Figure 5B). Next, to investigate specific cellular interactions in our dataset, we generated a list of the top ligand-receptor interactions for tumor/stroma cells versus immune cells ( Figure 5C), tumor/stroma cells versus tumor/ stroma cells (Figure S5A), and immune cells versus immune cells ( Figure S5B). Our analysis detected genes well studied in the NB field in addition to other interesting interactions. The genes included the presence of NOTCH pathway ligand and receptor genes known to play a role in NB tumor cell differentiation 72,73 and the expression of NRTK1 specifically in tumor cells, which, when activated, is described to cause neuronal differentiation and cell-cycle arrest in NB. 74 Subsequently three particularly interesting interactions (based on highest expression) drew our attention due to their relation to already described importance in cancer and correlation to survival. These may therefore be used as possible modulatory targets for NB treatment: semaphorin 6D (SEMA6D)/triggering receptor expressed on myeloid cells 2 (TREM2), galectin-9 (LGALS9)/ hepatitis A virus cellular receptor 2 (HAVCR2), and cluster of differentiation 24 (CD24)/sialic acid binding Ig-like lectin 10 (SI-GLEC10) ( Figure 5C). Survival analysis of the indicated interactions showed that high expression of SEMA6D-TREM2 was significantly correlated with decreased survival (Figure 5D), suggesting this interaction as a possible target in NB. On the contrary, high LGALS9-HAVCR2 expression was positively correlated with survival ( Figure 5E). In a similar manner, high CD24-SIGLEC10 expression was significantly correlated with improved survival ( Figure 5F). This suggests that LGALS9-HAVCR2 and CD24-SIGLEC10 might possibly be used in a different manner than the classical ''inhibition to unleash activity'' strategy. 75 Understanding and exploiting these interactions will require additional studies.

DISCUSSION
Here, we created a comprehensive cellular and molecular map of immune cells present within the NB microenvironment and provided a resource for future research into the development of novel therapeutics. We used a combination of single-cell transcriptomics, multiplex immunohistochemistry, and in silico analysis. The cell atlas shows 27 different immune cell subtypes that create a highly unique interactive network of cellular and molecular interactions within immune populations and with tumor and stroma cell populations in NB.
It is important to emphasize here that analyzing scRNA-seq from different studies where different tissue dissociation protocols were used might introduce bias in terms of cell fraction changes. We were therefore cautious when interpreting cell fraction differences since distortions from protocols are still likely to be pronounced. Instead, we focused our analyses on cell-state changes that should provide a more accurate analysis. 32 In addition, it is important to note that 7 out of the 17 patients included in our study were pre-treated and/or treated with systemic therapy at the time of tissue harvest. Before tissue harvest, patients had at least 2 weeks of wash out after chemotherapy and none of the patients were treated with immunotherapy. Recently, the first study trying to understand the influence of therapy on immune cells in the blood of high-risk NB patients was conducted. Interestingly, in this study, large variations between patients existed and immune cell phenotypes related to both tumor suppression and support were detected. 76 The influence of different cancer therapies on immune cells in the broader cancer perspective has not been defined yet, with chemotherapy having both immune enhancing and inhibiting properties. 77 Furthermore, interpretation of survival analysis can be difficult since the outcome is highly dependent on the type of analysis, the number of patient samples, the distribution of outcome, which biomarkers are analyzed, and possible additional factors. For several cell types, we detected a significant correlation with improved survival overall. However, when splitting the samples into different risk group and MYCN status the differences disappear. This may be due to the fact that low-risk patients do well overall; therefore, the number of high-risk patients goes down and it is difficult to detect differences. Other factors than tumor biology (like treatment) may influence outcome within the different subgroups (risk group/MYCN status); however, the limited number of patients did not allow us to study those. Nevertheless, the overall analysis and interpretation we have done is important for clinical purposes.
Myeloid cell infiltration has been described for multiple cancers and shown to support tumor growth. 18,41,42 In NB, CSF1R+ myeloid cells have been shown to predict poor clinical outcome 20 ; however, other studies have indicated myeloid cells to be correlated with improved survival. 78 The different myeloid cell subtypes in our study show a gene expression pattern indicating ongoing plasticity in both monocytes and macrophages, and we detected cells with pro-inflammatory phenotypes, which is in contrast to previous studies regarding myeloid cells in NB. 20,78 Pro-inflammatory immune cells and inflammation can promote tumor growth and cause metastatic spread in cancer. 79 However, our data suggest that, in NB, a pro-inflam-matory state is a good prognostic factor due to the positive correlation between different myeloid cell populations with pro-inflammatory cell states and improved survival. The different myeloid populations may still be exploited in the pursuit for novel therapeutics. Several interactions exist between tumor/stroma and myeloid cells, which could be exploited as therapeutic targets to recruit the good prognostic monocyte and macrophage players into the NB tumor. For this, we suggest further in-depth studies that will hopefully demonstrate their specific function and the possible use of these interactions in therapeutic strategies.
B cells have gained significantly more attention in the cancer field. Three studies demonstrated that B cells and the presence of TLSs in adult tumors correlated with improved response to immunotherapy in melanoma and sarcoma. 51,80,81 In human NB, a paper on the interaction between dendritic cells and NK cells detected the presence of TLSs. 14 Since TLSs are dynamic structures containing different immune cell types, they may be highly interesting for targeting or supporting therapy.
T cells, mainly CTLs, have been extensively studied in cancer. CTLs can play an important role in cancer eradication and have therefore been the main targets of currently used immunotherapy such as PD-L1 checkpoint inhibitors. 82 In our study, we showed that different CTLs are an important part of the interacting NB microenvironment and that several CTL subtypes correlated with improved survival. Importantly, we did not detect differences in exhaustion score across NB risk groups, implying limited phenotypic changes over the course of disease. Moreover and surprisingly, Tregs were not correlated with survival, a finding that is in contrast with Treg suppression and exhaustion of effector cells observed in adult cancer. 83 This indicates that the NB tumor microenvironment differs from adult cancer. Another possible explanation might be a result of dilution in bulk RNA sequencing from survival datasets to the low number of Tregs in the tumor microenvironment.
Our receptor-ligand interaction analysis between tumor/ stroma and immune cells revealed several genes that have already been correlated with NB growth and differentiation, [72][73][74]84 as well as the discovery of important cell-cell interactions mediating immune suppression in other cancer types, such as the CCL20-CCR6 axis in prostate cancer bone metastasis 32 and the CD161-CLEC2D axis in glioma. 85 This provides a basis for the additional genes detected in our analysis to have a potential therapeutic impact. For instance, SEMA6D was expressed mainly on SCP-like cells and functions as an inhibitory ligand upon binding to Plexin A1, which forms a receptor complex with TREM2, a marker for TAMs. 86 TREM2 deficiency has been shown to attenuate tumor growth, and combining TREM2 deficiency with PD-1 antibodies caused tumors to completely regress in a sarcoma mouse model. 87 Since we detected that SEMA6D-TREM2 expression is significantly correlated with decreased survival suggesting, we suggest this interaction as a potential therapeutic target. Furthermore, we detected elevated LGALS9 expression in endothelial cells that interacts with HAVCR2 (encodes TIM-3) on NK cells. TIM-3 present on NK cells serves as a maturation marker and functions as an inhibitory receptor, while TIM-3 induced inhibition can be overcome by 10 Cell Reports Medicine 3, 100657, June 21, 2022 Article ll OPEN ACCESS cytokine stimulation. 88 The function of TIM-3 in NK cells upon galectin-9 binding and immune cues differs depending on the cues and the endurance of the input, rendering TIM-3 function in NK cells in different cancers controversial. It has been suggested that the initial upregulation of TIM-3 enhances cytotoxicity and that chronic overexpression or dysregulation leads to exhaustion. 89 However, the precise role of TIM-3 in NK cell function in NB cannot be determined based on our analysis, and therefore warrants further investigation. Finally, the surface receptor CD24 was exclusively expressed on tumor cells and interacts with SIGLEC10 on macrophages in the immune cells. CD24 is a ''don't eat me'' signal, which conveys an inhibitory response in macrophages when binding to SIGLEC10 preventing these macrophages from exerting phagocytosis. The presence of CD24 has been detected on several solid tumors and its importance in suppressing phagocytosis has recently been demonstrated in ovarian and breast cancer. In addition, it is suggested that CD24/SIGLEC10 might be a target for cancer immunotherapy. 90 Hence, monoclonal anti-CD24 antibodies could be investigated as a potential therapy in NB.
Immune cells and immunotherapy in NB are an active field of research. Despite lack of MHC class I and high mutational load in NB tumors, major steps toward improved therapy are being made. Two recent studies demonstrated the potential of immunotherapy in NB by improving current anti-GD2 to reduce pain and by creating novel chimeric antigen receptors(CARs) on T cells toward a peptide sequence derived from the master transcriptional regulator in NB PHOX2B. 23,91 Our study provides a comprehensive resource for possible ways the tumor microenvironment could be exploited to further enhance immunotherapy in NB.
In conclusion, our study suggests that the immunological landscape in the NB microenvironment appears and most probably functions differently from the well-characterized adult cancer tumor microenvironment. Since NB is known to be initiated early during embryonal development, 1,31 it may be possible that tumor development occurred alongside immune cell development and therefore other mechanisms are at play within the NB microenvironment. We believe that our findings will pave the way to a better understanding of the role of the immune system in NB. It is important to emphasize that future studies are needed to understand the functional properties of each immune cell subtype and the suggested cell-cell interactions between the different cell types that our study has identified. Combined, we hope that these biological insights will help to tailor new and improved therapeutics that instruct the immune system against NB.
Limitations of the study While the presented analysis provides a good representation of immune cell heterogeneity in NB in general, a number of potential limitations should be noted. First, our dissociation protocol used for preparing single-cell suspensions is known to enrich for immune cells. 32 We, therefore, observe a higher frequency of immune cells relative to tumor and stromal cells in our dataset compared with public datasets and were cautious about interpreting frequency data outside of the data obtained using the same protocol. This enrichment, however, gave us the opportunity to study immune cell granularity in NB in great detail. Second, we have performed interpretation of functionality of certain immune cell types based on gene expression data. Protein validation and additional functional studies will have to provide additional evidence for these interpretations. Finally, creation and interpretation of survival curves based on gene expression can be challenging. NB is a highly heterogeneous disease in which genetics, risk group, age, and treatment can strongly influence outcome, and gene expression values are not the optimal choice for assessing survival. Despite the above-mentioned reasons, we managed to detect significant differences, which we believe provide important clinical insights into the influence of immune cells on NB survival.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Article ll Material availability Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ninib Baryawno (n.baryawno@ki.se).
Data and code availability d The expression datasets generated in this study are available through Gene Expression Omnibus with the accession number GEO: GSE147766. d Interactive views of the single cell datasets, differential expression results, code notebooks, cell annotation and RData objects are available on the author's website at https://github.com/shenglinmei/NB.immune.atlas/. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Patient material and collection of tumor specimens NB samples from Swedish patients were collected in conjunction with routine clinical sampling after parents and/or legal guardians had given verbal or written consent. All samples were collected according to permits approved by the Karolinska Institutet/Karolinska University Hospital ethics committees (reference numbers 03-736, 2009/1369-31/1) in accordance with the Helsinki declaration. Clinical management and therapy were performed according to national and international protocols and data were obtained from hospital records.
Tumor material used for immunohistochemistry 43 patients between the age 0 and 11 years were included in this study after written consent from their legal guardians. All patients have been treated at Lund University hospital from the year 2000 and forward. The tumor samples were accessed through the biobank at department of clinical genetics and pathology at Skå ne University hospital and regional pathology archives and consists of multiple, formalin-fixed paraffin embedded (FFPE) tumor samples.
Ethical approval for material used in immunohistochemistry All tissue samples and genetic data used during this project were collected from human tumors which were collected with informed consent from the legal guardians of the patients. Privacy of the donors have been protected and ethical approval from the ethical committee at Lund University has been obtained before the start of this project.

METHOD DETAILS
Dissociation of NB into single cells Fresh tissue obtained from surgery was collected in Media 199 supplemented with 2% (v/v) FBS on ice. Single cell suspensions of the tumors were obtained by cutting the tumor into small pieces (1mm 3 ) followed by enzymatic dissociation for 45 minutes at 37 C with shaking at 120 rpm using Collagenase I, Collagenase II, Collagenase III, Collagenase IV (all at a concentration of 1 mg/mL) and Dispase (2 mg/mL) in the presence of DNAse I solution (1:100) and RNase Out (1:1000). Cells were then resuspended in fetal bovine serum (FBS) with 5% DMSO and cryopreserved in liquid nitrogen.

FACS enrichment of viable NB cells
Cryopreserved cells from NB samples were thawed and stained with anti-CD235-PE and Calcein AM for 30 minutes at 4 C. Cells were washed twice with Media199 containing 2% FBS followed by DAPI and/or 7AAD staining (1 mg/mL). Flow sorting for live and non-erythroid cells (DAPI/7AAD-negative, Calcein AM-positive, CD235-negative) was performed on a BD FACS Aria III equipped with a 100mm nozzle instrument.

Single cell library preparation and sequencing
After FACS sorting, NB single-cell libraries were prepared using the Chromium single-cell 3 0 reagent kit v2 and v3 according to the manufacturer's recommendations. Libraries were sequenced on the Illumina NextSeq 500 platform.

Sectioning of NB tumors
Snap frozen NB samples were slowly brought from -80 C to +4 C. Samples were fixed in 4% paraformaldehyde in PBS (pH 7.4) at 4 C for 4 hours with constant rotation, rinsed in PBS 3 times and cryoprotected by incubating at 4 C overnight in PBS containing 30% sucrose, again with constant rotation. Tissue samples were then embedded in OCT and frozen at À20 C. Serial sections (10, 14 or 40 mm) were produced from each sample and collected on SuperFrost Plus Adhesion microscope slides and stored at À20 C. Highly multiplexed immunofluorescence For highly multiplexed immunofluorescence of tissue sections with the CODEX (Co-detection by inDEXing) we used the service provided by the Spatial Proteomics Facility at Scilifelab. We used the automated image acquisition with an inverted fluorescent microscope, Leica Dmi8, and fluidics exchange was performed using an Akoya CODEX instrument and CIM (CODEX Instrument Management software). Staining was done for the following targets: CD3, CD4, CD8, CD11c, CD45RO, CD19, HLA-DR and Ki67. A nuclear stain was used as a reference marker between the different cycles (SKU 7000003 CODEX Nuclear Stain).

Immunohistochemistry
Immunohistochemistry was performed on 4 mm tissue sections from FFPE tumor samples. The slides were dried for 20 min, then deparaffinized in Xylene and rehydrated in ethanol to water. Endogenous peroxidase was blocked for 20 min with 1% H2O2. The slides were then washed, and epitope retrieval was performed using Target Retrieval Solution (DAKO) in 95 C in a pressure cooker heater for 20 min. Slides were blocked in Protein Block Serum-Free Solution (X0909, DAKO) and later incubated with a primary antibody diluted in PBS containing 5% normal goat serum. Following incubation with primary antibody slides were washed in PBS and incubated with a secondary antibody. Labeled-chromogen substrate solution with Liquid DAB chromogen was applied. The slides were counterstained in Mayer htx before rehydrated and mounted with Faramount Aqueous Mounting Medium (S3025, DAKO).

Infiltration density of immune cells in human NB
The infiltration density of immune cells in the tumors were evaluated by counting cells in 4 mm tissue sections from FFPE tumor samples. The FFPE samples had been stained using immunohistochemistry for the cell markers; CD68, CD3, CD8, FOXP3, IL17, NE, NKp46 and CD1a. 100 cells were counted in three different areas in each sample at 40x magnification (Olympus BX43 microscope, Olympus, Shinjuku, Japan) and the average percentage of infiltrating cells were calculated for each sample.

QUANTIFICATION AND STATISTICAL ANALYSIS
Pre-processing of scRNA-seq data Demultiplexing of bcl files into fastq files was performed using cellranger 3.0.1 mkfastq software and alignments to the human genome reference sequences were performed using cellranger count. The reference included all protein coding genes as well as mitochondrial genes for downstream analysis. All cell barcodes with less than 500 UMIs were excluded. These were further filtered one sample at a time where barcodes with percent mitochondrial reads larger than median plus two standard deviations (percent.mito > median + 2sd). Likewise, barcodes with few detected genes were filtered out as log 10 (nGene) < median -2sds.

ScRNA-seq analysis
During the preliminary round of analyses the datasets were aligned using Conos 1.2.1 with k = 15; kself = 5 PCA rotation space and angular distance measure. Visualization using UMAP embedding showed a number of continuous bridges connecting the major populations.
For subtype assessment within myeloid, T, NK or B cells, we extracted all myeloid (or T, NK, B cells), removed low quality samples with less than 50 cells per cells, and realigned separately using Conos with default paramenters. Leiden community clustering method (as implemented in Conos) was used to determine refined joint clusters, providing higher resolution than the initial analysis. In addition, we also provided walktrap.community and multilevel.community result in data website (https://github.com/shenglinmei/ NB.imzmune.atlas), where users can download and view different clustering results. To avoid over clustering, we also evaluated cell type or cell states specific expressed genes. Final cell annotations are annotated based on literature reported markers and cluster specific markers (Table S5). getDifferentialGenes function from Conos was used to identify differentially expressed (marker) genes for clusters or subtypes. Annotation of the cluster communities was done using marker gene expression. Two patient samples were excluded for immune cell analysis due to low number of immune cells. These samples were included in the receptor ligand analysis where tumor and stroma are taken into account (see section Identification of significant ligand-receptor pairs).

Analysis of cell composition shift
We present two methods to analysis relative cell compositions shift within the major cell populations. The first one is clusterbased cell proportion analysis, cell frequency is normalized by total number of cells per sample (like total myeloid cells), statistical significance of proportion differences was then evaluated using Wilcoxon rank sum test. The second one is clusterfree cell composition analysis, where we estimate cell density per sample and evaluate the differential cell density between sample groups. We first compute kernel density in joint embedding space for each sample using ks R package (bin = 400), then quantile normalization was used to normalize the density matrix across samples and the average density of each sample group is shown in Figure 2H Calculation of gene set signature scores To assess cell states in different cell subsets and conditions, we used a gene set signature score to measure the relative difference of cell states. The signature scores were calculated as average expression values of the genes in a given set. Specifically, we first calculated signature score for each cell as an average normalized (for cell size) gene expression magnitudes, then the signature score for each sample was computed as the mean across all cells.

Survival analysis
We collected 498 bulk RNA sequencing data of primary NB patient samples from GSE49711 (https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE49711). To test if a given signature predicts survival of NB patient, we first computed the average gene expression level of the signature in each tumor based on the bulk RNA sequencing data. The samples were grouped into high (top 25%) and low (bottom 25%) groups based on the average signature gene expression. Next, we used a one-sided log-rank test to compute the significance of the association between the signature expression value and prognosis. To evaluate the satiability of signature genes list, we resample the signature genes and repeat the analysis with 500 bootstrap resampling rounds. Statistical significance was then assessed by p-values at the 99% reproducibility power (i.e. reporting 0.99th quantile of the sampled p-values). Similarly, the same survival analysis was separately applied to low-risk and high-risk patients, MYCN amplified and MYCN nonamplified patients. The detailed genes list used for prognostic analysis can be found in Table S5. For survival heatmap plots, P-value from the survival plot was taken and presented as heatmap using -log10 P-value.

Identification of significant ligand-receptor pairs
The ligand-receptor (LR) pairs were downloaded from published databases CellPhoneDB (v1.1.0). To identify significant ligand-receptor pairs in 10X data (n = 19 patient samples), we used a similar approach as previously described (Vento-Tormo et al., 2018). We first calculated gene expression ratio in each cell type and only considered genes with more than 10% of cells demonstrating expression within each cell type. We then calculated average expression of ligand and receptor pairs across cell type pairs in normalized scRNA-seq data. The product of average ligand expression in cell type A and the average receptor expression in cell type B was used to measure LR pair expression. Statistical significance was assessed by randomly shuffling the cluster labels of all cell types and recalculating ligand-receptor average pair expression across 1,000 permutations, which generated a null distribution for each LR pair in each pairwise comparison between two cell types. P-values were calculated with the normal distribution curve generated from the permuted LR pair interaction scores. Candidate LR pairs were determined by P value 0.05 or lower. To prioritize the significant ligandreceptor pairs, we tested if a LR pair is associated with patient overall survival. We first computed the average gene expression of the ligand and receptor in each tumor based on the bulk RNA-Seq data, next patients were stratified into two groups according to the average expression of the RL pair: high or low expression correspond to the top or bottom 25% of the population, respectively. Logrank test was used to examine if there was a significant difference between these two patient groups in terms of their survival. LR pairs were filtered using survival significance p value of 0.05 or lower. In addition, we also evaluated ligand and receptor expression, requiring both ligand and receptor highly expressed in corresponding cell type. The getDifferentialGenes() function from Conos was used to derive DEG from each cell type and genes. We next screened each of the LR pair using p-value determined Z score >5 and AUC > 0.5. The detailed ligand-receptor pairs can be found at: https://github.com/shenglinmei/NB.immune.atlas.

Statistics
Wilcoxon rank sum test was used to assess significance in signature score analysis and cell proportion differences analysis. The p-values in the figures were reported using the following symbols: *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001.

Data and code availability
Raw sequencing and processed data in this paper are available under the accession number: GSE147766. For the joint alignment analysis with public scRNA-seq data, we downloaded normal fetal adrenal tissue and NB scRNA-seq datasets from Dong et al., 2020, Kameneva et al., 2021 and Kildisiute et al., 2021. Custom code and the combined data that was used in this study can be found at https://github.com/shenglinmei/NB.immune.atlas. In addition, we provided an interactive view Conos object, allowing the user to download and view the data to the Conos joint alignment results.